/*
 * Copyright (c) 2022 Huawei Device Co., Ltd.
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

class NavierStokes {
    public static FluidField solver;
    public static int nsFrameCounter = 0;
    public static double factor = 0.5;
    public static int densFactor = 10;
    public static int point = 64;
    public static int repeatCount = 15;
    public static int checkIndex = 7000;
    public static int checkIndexTo = 7100;
    public static int indexTwo = 2;
    public static int point2 = 128;
    public static int resultOK = 77;
    public static int densInt0 = 5;
    public static int densInt1 = 20;
    public static int densInt2 = 30;
    public static int itersMax = 100;
    public static int pointMax = 1000000;
    public static int solveNumber = 4;
    public static int multiple = 1000000;
    public static int runCount = 120;
    
    public static void runNavierStokes() {
        solver.update();
        nsFrameCounter += 1;
        if (nsFrameCounter == repeatCount) {
            checkResult(solver.getDens());
        }
    }

    public static void checkResult(double[] dens) {
        int result = 0;
        for (int i = checkIndex; i < checkIndexTo; i++) {
            result += (int) (dens[i] * densFactor);
            // Log.print("single result[" + i + "] = " + (int) (dens[i] * densFactor));
        }
        // Log.print("checkResult result = " + result);
        if (result != resultOK) {
            // Log.print("nsFrameCounter = " + nsFrameCounter + " checksum failed");
        } else {
            // Log.print("nsFrameCounter = " + nsFrameCounter + " checksum success");
        }
    }
    /*
     * @Setup
     */
    public static void setupNavierStokes() {
        solver = new FluidField();
        solver.setResolution(point2, point2);
        solver.setIterations(densInt1);
        solver.reset();
    }
    public static int framesTillAddingPoints = 0;
    public static int framesBetweenAddingPoints = 5;

    public static void prepareFrame(double[] d, double[] u, double[] v) {
        if (framesTillAddingPoints == 0) {
            solver.addPoints(d, u, v);
            framesTillAddingPoints = framesBetweenAddingPoints;
            framesBetweenAddingPoints += 1;
        } else {
            framesTillAddingPoints -= 1;
        }
    }
}

class FluidField {
    public void addFields(double[] x, double[] s, double dt) {
        for (int i = 0; i < size; i++) {
            x[i] += dt * s[i];
        }
    }

    public void setBnd(int b, double[] x) {
        if (b == 1) {
            for (int i = 1; i <= width; i++) {
                x[i] = x[i + rowSize];
                x[i + (height + 1) * rowSize] = x[i + height * rowSize];
            }
			
            for (int j = 1; j <= height; j++) {
                x[j * rowSize] = -x[1 + j * rowSize];
                x[width + 1 + j * rowSize] = -x[width + j * rowSize];
            }
        } else if (b == NavierStokes.indexTwo) {
            for (int i = 1; i <= width; i++) {
                x[i] = -x[i + rowSize];
                x[i + (height + 1) * rowSize] = -x[i + height * rowSize];
            }
			
            for (int j = 1; j <= height; j++) {
                x[j * rowSize] = x[1 + j * rowSize];
                x[width + 1 + j * rowSize] = x[width + j * rowSize];
            }
        } else {
            for (int i = 1; i <= width; i++) {
                x[i] = x[i + rowSize];
                x[i + (height + 1) * rowSize] = x[i + height * rowSize];
            }
			
            for (int j = 1; j <= height; j++) {
                x[j * rowSize] = x[1 + j * rowSize];
                x[width + 1 + j * rowSize] = x[width + j * rowSize];
            }
        }
        int maxEdge = (height + 1) * rowSize;
        x[0] = NavierStokes.factor * (x[1] + x[rowSize]);
        x[maxEdge] = NavierStokes.factor * (x[1 + maxEdge] + x[height * rowSize]);
        x[width + 1] = NavierStokes.factor * (x[width] + x[width + 1 + rowSize]);
        x[width + 1 + maxEdge] = NavierStokes.factor * (x[width + maxEdge] + x[width + 1 + height * rowSize]);
    }

    public void linSolve(int b, double[] x, double[] x0, int a, int c) {
        if (a == 0 && c == 1) {
            for (int j = 1; j <= height; j++) {
                int currentRow = j * rowSize;
                currentRow += 1;
                for (int i = 0; i < width; i++) {
                    x[currentRow] = x0[currentRow];
                    currentRow += 1;
                }
            }
            setBnd(b, x);
        } else {
            double invC = 1.0 / c;
            for (int k = 0; k < iterations; k++) {
                linSolveElse(x, x0, a, invC);
                setBnd(b, x);
            }
        }
    }
    public void linSolveElse(double[] x, double[] x0, int a, double invC) {
        for (int j = 1; j <= height; j++) {
            int lastRow = (j - 1) * rowSize;
            int currentRow = j * rowSize;
            int nextRow = (j + 1) * rowSize;
            double lastX = x[currentRow];
            int changeRow = currentRow;
            currentRow += 1;
            for (int i = 1; i <= width; i++) {
                double args1 = x0[currentRow];
                changeRow = currentRow;
                currentRow += 1;
                lastRow += 1;
                nextRow += 1;
                double args2 = a * (lastX + x[currentRow] + x[lastRow] + x[nextRow]);
                x[changeRow] = (args1 + args2) * invC;
                lastX = x[changeRow];
            }
        }
    }

    public void diffuse(int b, double[] x, double[] x0) {
        int a = 0;
        linSolve(b, x, x0, a, 1 + NavierStokes.solveNumber * a);
    }

    public void linSolve2(double[] x, double[] x0, double[] y, double[] y0, int a, int c) {
        if (a == 0 && c == 1) {
            for (int j = 1; j <= height; j++) {
                int currentRow = j * rowSize;
                currentRow += 1;
                for (int i = 0; i < width; i++) {
                    x[currentRow] = x0[currentRow];
                    y[currentRow] = y0[currentRow];
                    currentRow += 1;
                }
            }
            setBnd(1, x);
            setBnd(NavierStokes.indexTwo, y);
        } else {
            double invC = 1.0 / c;
            for (int k = 0; k < iterations; k++) {
                linSolve2Else(x, x0, y, y0, a, invC);
                setBnd(1, x);
                setBnd(NavierStokes.indexTwo, y);
            }
        }
    }

    public void linSolve2Else(double[] x, double[] x0, double[] y, double[] y0, int a, double invC) {
        for (int j = 1; j <= height; j++) {
            int lastRow = (j - 1) * rowSize;
            int currentRow = j * rowSize;
            int nextRow = (j + 1) * rowSize;
            double lastX = x[currentRow];
            double lastY = y[currentRow];
            currentRow += 1;
            for (int i = 1; i <= width; i++) {
                currentRow += 1;
                lastRow += 1;
                nextRow += 1;
                double test = x0[currentRow] + a * (lastX + x[currentRow] + x[lastRow] + x[nextRow]);
                x[currentRow] = test * invC;
                lastX = x[currentRow];
                double testy = y0[currentRow] + a * (lastY + y[currentRow] + y[lastRow] + y[nextRow]);
                y[currentRow] = testy * invC;
                lastY = y[currentRow];
            }
        }
    }

    public void diffuse2(double[] x, double[] x0, double[] y, double[] y0) {
        int a = 0;
        linSolve2(x, x0, y, y0, a, 1 + NavierStokes.solveNumber * a);
    }

    public void advect(int b, double[] d, double[] d0, double[] u, double[] v, double dt) {
        double wdt0 = dt * width;
        double hdt0 = dt * height;
        double wp5 = width + NavierStokes.factor;
        double hp5 = height + NavierStokes.factor;
        for (int j = 1; j <= height; j++) {
            int pos = j * rowSize;
            for (int i = 1; i <= width; i++) {
                pos += 1;
                double x = i - wdt0 * u[pos];
                double y = j - hdt0 * v[pos];
                if (x < NavierStokes.factor) {
                    x = NavierStokes.factor;
                } else if (x > wp5) {
                    x = wp5;
                }
                int i0 = (int) x | 0;
                int i1 = i0 + 1;
                if (y < NavierStokes.factor) {
                    y = NavierStokes.factor;
                } else if (y > hp5) {
                    y = hp5;
                }
                int j0 = (int) y | 0;
                int j1 = j0 + 1;
                double s1 = x - i0;
                double s0 = 1 - s1;
                double t1 = y - j0;
                double t0 = 1 - t1;
                int row1 = j0 * rowSize;
                int row2 = j1 * rowSize;
                d[pos] = s0 * (t0 * d0[i0 + row1] + t1 * d0[i0 + row2])
                        + s1 * (t0 * d0[i1 + row1] + t1 * d0[i1 + row2]);
            }
        }
        setBnd(b, d);
    }

    public void project(double[] u, double[] v, double[] p, double[] div) {
        double h = -NavierStokes.factor / Math.sqrt(width * height);
        for (int j = 1; j <= height; j++) {
            int row = j * rowSize;
            int previousRow = (j - 1) * rowSize;
            int prevValue = row - 1;
            int currentRow = row;
            int nextValue = row + 1;
            int nextRow = (j + 1) * rowSize;
            for (int i = 1; i <= width; i++) {
                currentRow += 1;
                nextValue += 1;
                nextRow += 1;
                previousRow += 1;
                prevValue += 1;
                div[currentRow] = h * (u[nextValue] - u[prevValue] + v[nextRow] - v[previousRow]);
                p[currentRow] = 0;
            }
        }
        setBnd(0, div);
        setBnd(0, p);
        linSolve(0, p, div, 1, NavierStokes.solveNumber);
        double wScale = NavierStokes.factor * width;
        double hScale = NavierStokes.factor * height;
        for (int j = 1; j <= height; j++) {
            int prevPos = j * rowSize - 1;
            int currentPos = j * rowSize;
            int nextPos = j * rowSize + 1;
            int prevRow = (j - 1) * rowSize;
            int nextRow = (j + 1) * rowSize;
			
            for (int i = 1; i <= width; i++) {
                currentPos += 1;
                nextPos += 1;
                prevPos += 1;
                nextRow += 1;
                prevRow += 1;
                u[currentPos] -= wScale * (p[nextPos] - p[prevPos]);
                v[currentPos] -= hScale * (p[nextRow] - p[prevRow]);
            }
        }
        setBnd(1, u);
        setBnd(NavierStokes.indexTwo, v);
    }

    public void densStep(double[] x, double[] x0, double[] u, double[] v, double dt) {
        addFields(x, x0, dt);
        diffuse(0, x0, x);
        advect(0, x, x0, u, v, dt);
    }

    public void velStep(double[] u, double[] v, double[] u0, double[] v0, double dt) {
        addFields(u, u0, dt);
        addFields(v, v0, dt);
        double[] temp = u0;
        u0 = u;
        u = temp;
        double[] temp1 = v0;
        v0 = v;
        v = temp1;
        diffuse2(u, u0, v, v0);
        project(u, v, u0, v0);
        double[] temp2 = u0;
        u0 = u;
        u = temp2;
        double[] temp3 = v0;
        v0 = v;
        v = temp3;
        advect(1, u, u0, u0, v0, dt);
        advect(NavierStokes.indexTwo, v, v0, u0, v0, dt);
        project(u, v, u0, v0);
    }

    public void queryUI(double[] d, double[] u, double[] v) {
        for (int i = 0; i < size; i++) {
            u[i] = 0.0;
            v[i] = 0.0;
            d[i] = 0.0;
        }
        NavierStokes.prepareFrame(d, u, v);
    }

    public void update() {
        queryUI(densPrev, uPrev, vPrev);
        velStep(u, v, uPrev, vPrev, dt);
        densStep(dens, densPrev, u, v, dt);
    }

    public void setIterations(int iters) {
        if (iters > 0 && iters <= NavierStokes.itersMax) {
            iterations = iters;
        }
    }

    FluidField() {
        iterations = 10;
        visc = 0.5;
        dt = 0.1;
        setResolution(NavierStokes.point, NavierStokes.point);
    }

    public int iterations;
    public double visc;
    public double dt;
    public double[] dens = new double[0];
    public double[] densPrev = new double[0];
    public double[] u = new double[0];
    public double[] uPrev = new double[0];
    public double[] v = new double[0];
    public double[] vPrev = new double[0];
    public int width = 0;
    public int height = 0;
    public int rowSize = 0;
    public int size = 0;
    public void reset() {
        rowSize = width + NavierStokes.indexTwo;
        size = (width + NavierStokes.indexTwo) * (height + NavierStokes.indexTwo);
        dens = new double[size];
        densPrev = new double[size];
        u = new double[size];
        uPrev = new double[size];
        v = new double[size];
        vPrev = new double[size];
        for (int i = 0; i < size; i++) {
            densPrev[i] = uPrev[i] = vPrev[i] = dens[i] = u[i] = v[i] = 0;
        }
    }

    public void setResolution(int hRes, int wRes) {
        int res = wRes * hRes;
        if (res > 0 && res < NavierStokes.pointMax && (wRes != width || hRes != height)) {
            width = wRes;
            height = hRes;
            reset();
        }
    }

    public double[] getDens() {
        return dens;
    }

    public void addPoints(double[] md, double[] mu, double[] mv) {
        int n = 64;
        for (int i = 1; i <= n; i++) {
            setVelocity(mu, mv, i, i, n, n);
            setDensity(md, i, i, NavierStokes.densInt0);
            setVelocity(mu, mv, i, n - i, -n, -n);
            setDensity(md, i, n - i, NavierStokes.densInt1);
            setVelocity(mu, mv, NavierStokes.point2 - i, n + i, -n, -n);
            setDensity(md, NavierStokes.point2 - i, n + i, NavierStokes.densInt2);
        }
    }

    public void setDensity(double[] den, int x, int y, double d) {
        den[(x + 1) + (y + 1) * rowSize] = d;
    }

    public void setVelocity(double[] u, double[] v, int x, int y, int xv, int yv) {
        u[(x + 1) + (y + 1) * rowSize] = xv;
        v[(x + 1) + (y + 1) * rowSize] = yv;
    }
}

/*
 * @State
 * 
 * @Tags Jetstream2
 */
class Benchmark {

    public static void main(String[] args) {
        new Benchmark().run();
    }

    Benchmark() {
        NavierStokes.setupNavierStokes();
    }

    public void runIteration() {
        NavierStokes.runNavierStokes();
    }
    /*
     * @Benchmark
     */
    public void run() {
        long start = System.nanoTime();
        for (int i = 0; i < NavierStokes.runCount; i++) {
            runIteration();
        }
        long end = System.nanoTime();
        System.out.println("navier-stokes: ms = " + (double)(end - start) / NavierStokes.multiple);
    }
}

class Log {
    static boolean debug = false;

   public static void print(String msg) {
        if (debug) {
            System.out.println(msg);
        }
    }

}